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Abstract 

We report the results of a detailed numerical study designed to estimate both 
the absolute age and the uncertainty in age (with confidence limits) of the 
oldest globular clusters. Such an estimate is essential if a comparison with 



c§ ! the Hubble age of the universe is to be made to determine the consistency, 

or lack thereof, of various cosmological models. Utilizing estimates of the 
^ \ uncertainty range (and distribution) in the input parameters of stellar evolu- 

tion codes we produced 1000 Monte Carlo realizations of stellar isochrones, 
with which we could fit the ages of the 18 oldest globular clusters. Incor- 
porating the observational uncertainties in the measured color-magnitude 
diagrams for these systems and the predicted isochrones, we derived a prob- 
ability distribution for the mean age of these systems. The one-sided 95% 
C.L. lower bound for this distribution occurs at an age of 12.07 Gyr. This 
puts interesting constraints on cosmology which we discuss. Further details, 
including a description of the distributions, covariance matrices, dependence 
upon individual input parameters, etc. will appear in a future article. 



^Also Dept of Astronomy. 
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The apparent dichotomy between the upper bound on the age of the 
Universe obtained from the Hubble Constant and the lower bound obtained 
by dating the oldest globular clusters in our galactic halo represents one of 
the most significant potential conflicts in modern observational cosmology. 
For Hubble constant Hq = 100 h kmsec~^Mpc ^, a flat, matter dominated 
universe has an age given by 

THubbie = 2/3 H^^ = Q.Qh-'Gyr (1) 

If the universe is open, the factor of 2/3 is changed to unity, but even in 
this case if the age of the oldest globular clusters in our galaxy is really 16 
Gyr, as various best fits estimates recently suggest |l|, ||, then this will be 
inconsistent with the Hubble age ii h > .6. Consistency with the fiat matter 
dominated universe estimate is impossible ii h > A. Since several recent 
estimates of the Hubble constant based on Type la supernova, and Hubble 
Space Telescope observations of Cepheid variables in the Virgo cluster both 
tend to suggest /i > .65 this has been one factor which has led various 

groups to argue once again for the need for a cosmological constant (i.e. |^). 

Since the Hubble estimate is unambiguous for a fixed Hubble constant the 
crucial uncertainty in this comparison resides in the globular cluster (GC) 
ages estimates themselves. Rough arguments have been made that changes 
in various input parameters in the stellar evolution codes designed to derive 
globular cluster isochrones, or in the RR Lyrae distance estimator used to 
determine absolute magnitudes for GC stars, might change age estimates by 
10-20 % (i.e. p). However, no systematic study has yet been undertaken to 
realistically estimate the cumulative effect of all existing observational and 
theoretical uncertainties in the GC age analysis. This is the purpose of the 
present work. 
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One of the reasons we believe such an analysis had not yet been carried 
out is that it is numerically intensive. Each run of a stellar evolution code for 
a single mass point takes 3-5 minutes on the fastest commercially available 
workstations. Nine different mass points at three different metallicity values 
must be run to produce each set of isochrones. If one then runs, say, 1000 
different isochrone sets to explore the different parameter ranges available, 
this requires over 8 weeks of continuous processing time. 

This is long, but not prohibitive, so because of the importance of this 
issue, we developed the Monte Carlo algorithms necessary for the task. This 
involved first examining the measurements of input parameters in the stellar 
evolution code to determine their best fit values, and also their uncertainties 
along with the appropriate distributions to use in the Monte Carlo. Then 
the stellar evolution code 0] and isochrone generation code were rewritten 
to allow sequential input of parameters chosen from these distributions, and 
output of the necessary color-magnitude (CM) diagram observables. Finally, 
we derived a fitting program to compare the predictions to the data. Since 
the numerically intensive part of this procedure involves the Monte Carlo 
generation of isochrones, by incorporating the chief observational CM un- 
certainty afterwards our results can quickly be refined as this uncertainty is 
refined. 

In this article we briefly describe the general features of our analysis and 
present our main result, the derived probability distribution for the age of 
the oldest galactic globular clusters. This should be the result which is of 
broadest and most immediate general interest. In a future work we shall 
describe the details of the analysis, and present our results for covariance 
matrices, correlation functions, dependence on individual parameters, and 
observational uncertainties. These results will be of more use to researchers 
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in the field of stellar evolution and globular cluster dating. 

Monte Carlo Analysis Inputs: General Features 

We have focussed on what we believe are the chief input uncertainties in 
the derivation of stellar evolution isochrones. These include: pp and CNO 
chain nuclear reaction rates, stellar opacity uncertainties, uncertainties in 
the treatment of convection and diffusion, helium abundance uncertainties, 
and uncertainties in the abundance of the a-capture elements (O, Mg, Si, 
S, and Ca). Our stellar evolution code was revised to allow batch running 
with sequential input of these parameters chosen from underlying probability 
distributions. We did not include the equation of state among our Monte 
Carlo variables as it is now well understood in metal-poor main sequence 
stars. It has recently been shown that the detailed equation of state by 
Rogers gives very similar globular cluster age estimates to those obtained 
using the Debye-Hiickel correction 0, which was used in this study. 

We describe briefly our parameter choices and distributions below: We 
included uncertainties for the 3 most important reactions in the pp chain. For 
p + p —^'^B.+e'^ + z/ we used the analysis of Kamionkowski and Bahcall (||1U||), 
except that where they took theoretical errors as gaussianly distributed, we 
decided that a uniform distribution better represented the state of our (lack 
of) knowledge. There are two sources of theoretical error, one from the 
uncertainty in the particle wave-functions, another from meson exchange 
111). We use a relative modification to this reaction of 1 zt .002 i'ooog J3i2^j 



where the 2nd term is 1-sigma gaussian, and the 3rd and 4th are top hat 
distributions. For 2 other pp chain reactions, ^He + ^He ^He + p + p, and 
^He + ^He ''Be + 7, the uncertainties included in this study were taken 
from Bahcall and Pinsonneault ([0, Table I). The CNO nuclear reaction 
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rates and their uncertainties are from Bahcall (||T3||; Table 3.4). 

The stellar opacities are broken into high temperature and low temper- 
ature regimes. For the high temperature regime (T > 10^ K), the OPAL 
opacities are utilized. The uncertainties in these opacities were evalu- 
ated by a comparison to the LAOL opacities [jl5|. In the temperature regime 
relevant for nuclear fusion {T ^6 x 10^ K), typically differences of 1% were 
found, with maximum differences of 3%. Thus, we elected to multiply the 
high temperatures opacities by a Gaussian distribution, with a mean of 1, 



and a = 0.01. The Kurucz opacities were used for the low-temperature 
regime. Low temperature opacities calculated by different groups can differ 



by a large amount fT^ , though modern calculations appear to agree to within 
~ 30%. In addition, we intercompared Kurucz's calculations for different el- 
ement mixtures, finding maximum differences of 30%. For these reasons, we 
elected to multiply the Kurucz opacities by a number which was uniformly 
drawn from the range 0.7 - 1.3. 

The correct treatment of the convective regions of the star has been a long 
standing problem in stellar astrophysics, which stems from our poor under- 
standing of highly compressible convection and its interaction with radiation 
in the optically thin outer layers. Most stellar evolution codes use the mix- 
ing length approximation, which is a simple model of convection where blobs 
of matter are supposed to rise or fall adiabatically over some distance, and 
then instantaneously release their heat |]18[. The uncertainties in this treat- 
ment of convection are parameterized into a single non-dimensional quantity 
called the mixing length, which is usually taken to be fixed during a star's 
evolution. Its value is typically chosen by requiring that a solar model have 
the correct radius and luminosity at the solar age, or by a comparison of 
theoretical isochrones to observed CM diagrams. A mixing length of around 
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1.8 appears to provide a reasonable match to the observations, but its ex- 
act value depends on the input physics (opacities, model atmospheres, etd.) 
used to construct the stellar models. Modern solar models typically employ 
mixing lengths between 1.7 and 2.1 |[19| , To further explore this issue, we 
have conducted tests whereby isochrones were constructed from models with 
different mixing lengths and compared to a metal-poor globular cluster CM 
diagram. We found that changes in the mixing length of 0.3 could be ruled 
out as the isochrones no longer fit the data, if all other input parameters 
were held constant. Allowing for possible variation among globular clusters 
and for the fact that a broader range in mixing length might fit the data if 
other parameters are allowed to vary at the same time, we elected to use a 
rather broad Gaussian distribution with a mean of 1.85 and a = 0.25. 

Whether or not to include the effects of element diffusion (whereby helium 
tends to sink to the center of the star, while hydrogen is rises to the surface) 
is a difficult question. Physical models of compressible plasmas suggest that 
diffusion should be occurring in stars, and predict diffusion coefficients with 
a claimed accuracy of about 30% |PT|. However, these models assume that 
no other mixing process occurs within the radiative regions of a star, which 
may not be true. Helioseismology appears to suggest that diffusion is occur- 



ring in the sun |2^, but the evidence is not compelling [|T9]. Models of halo 
stars which incorporate diffusion predict curvature in the Li-effective temper- 
ature plane which is not observed p3|, which suggests that some process is 
inhibiting diffusion in halo stars. Given these uncertainties, we have elected 



to multiply the diffusion coefficients given by Michaud & Proffitt [21] by a 
number uniformly drawn from the interval 0.3 - 1.2. The use of a flat distri- 
bution with a rather large range reflects our conviction that the incorporation 
of diffusion within stellar models is subject to large uncertainties. 
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The primordial helium abundance, relevant to old halo stars, and an 
important input in our stellar codes, has taken on renewed interest as a 
result of recent calculations of Big Bang Nucleosynthesis (BBN) light element 
production ^ In order to BBN estimates to agree with inferred 
primordial abundances, significant systematic uncertainties must be allowed 
for. It has become clear that such uncertainties are the dominant feature of 
the comparison between theory and observation. We therefore utilize here 
a flat distribution for the primordial helium mass fraction between 0.22 and 
0.25, which encompasses the range of recent estimates. 

In order to calculate a stellar model, the abundance of the elements heav- 
ier than helium (denoted by Z) must be specified. Due to its numerous 
spectral lines, it is relatively easy to determine the abundance of Fe in globu- 
lar cluster stars. Unfortunately, the abundances of the other heavy elements 
are more difficult to determine and it has been common to assume that the 
other heavy elements are present in the same proportion as they are in the 
Sun. However, from both theoretical arguments and observational evidence, 
it is clear that the elements which are produced via a-capture are enhanced 
in abundance relative to their solar value. It is relatively easy to incorporate 
the effects of the enhancement of the a-elements on the stellar models by re- 
defining the relationship between the model Z and the iron abundance [^]. 
Oxygen is by far the most important of the a-capture elements, as it acounts 
for roughly half of all the heavy elements (by number) present in the Sun. 
For this reason, we concentrate on observations of oxygen as representative 
of the a-capture elements. The determination of oxygen abundances in stars 
is extremely difficult, and subject to a number of systematic uncertainties 
(e.g. 10). Recently, high quality [O/Fe]0 abundances for a number of halo 
^We use the common spectroscopic notation, where the abundance of element y relative 
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stars have been obtained with the result that the mean abundance was 
found to be [0/Fe] = 0.55 ± 0.05, where the error is simply the standard 
deviations of the measurements. In addition to this error, one must add in 
the possible systematic errors which have been discussed by a number of 
authors An analysis of the literature has lead us to conclude 



that possible systematic errors in the determination of [0/Fe] may be as 
large as ±0.2 dex. Thus, the abundance of the a-elements was taken to be 
[a/Fe] = 0.55 ± 0.05± (Gaussian) ±0.2 (top-hat). 

Finally, a color table must be used to convert our theoretical luminosities 
and temperatures to observed magnitudes and colors. The construction of an 
accurate color table requires the use of theoretical model atmospheres, which 
are still subject to large uncertainties. We elected to take this uncertainty 
into account by randomly choosing one of two totally independent color ta- 
bles p2| , |33| with equal probability in constructing each isochrone set in the 



Monte Carlo. These two tables span reasonably the present range used to 
transform from theoretical temperatures and luminosities to observed colors 
and magnitudes. 

The Fitting Procedure and the Probability Distribution for Glob- 
ular Cluster Ages 

Once a set of isochrones (which consists of isochrones of three different 
metallicities for a set of 15 different ages between 8 and 22 Gyr) is derived, the 
comparison with the observed parameters of a specific set of globular clusters 
requires a fitting procedure. In order to minimize the large uncertainties in 



the effective temperatures of the models |34], we have elected to use the dif 



ference in magnitude between the main sequence turn-off, and the horizontal 



to element x is denote by [y/x] = log(?//a;)star - ^og{y/x)s 
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branch (HB, in the RR Lyr instabihty strip) as our age diagnostic. This 
age determination technique is commonly referred to as AVj'Jb ^-nd has been 
extensively used in the astronomical literature (e.g. Our isochrone sets 



provide us with the main-sequence turn-off luminosity as a function of age 
and metallicity. Due to the importance of convection in the nuclear burning 
regions of HB stars, theoretical HB luminosities are subject to large uncer- 
tainties, and so we have elected to combine our theoretical main-sequence 
turn-off luminosities with an observed relation for the luminosity of the HB 
(see below). This results in a grid of predicted AV^b's as a function of age 
and [Fe/H], which is then fit to an equation of the form 

t, = f3o + AAV + /?2AV2 + /?3[Fe/H] + f3,[Fe/}if + /?5AV[Fe/H], (2) 

where tg is the age in Gyr. The observed values of AV^b and [Fe/H], along 
with their corresponding errors, are input in (|^) to determine the age and its 
error for each GC in our sample. 

There is abundant evidence for a large age range within different GC 
systems (e.g. ||36|, ^) so one must take care to select a sample which only 



includes old globular clusters. Observational errors in the determination of 
the turn-off and horizontal branch magnitudes lead to a ~ 10 — 20% error 
in the derived age of any single cluster. Thus, to minimize the observa- 
tional uncertainties, it is best to determine the mean age of a number of 
GCs. In light of the strong evidence for an age-metallicity relationship (with 
metal-poor clusters being the oldest), only metal-poor clusters were selected 
([Fe/H] < —1.6). From this list of metal-poor clusters, any cluster which 
has been shown to be young using the difference in color between the giant- 
branch and main sequence turn-off [37], or which is suspected of being young 



due to its unusually red horizontal branch for its metallicity [^ was dis- 

9 



carded. From the sample of 43 GCs for which high quahty observations 
are available, 27 survived the metallicity cut, of which 10 were discarded for 
being young, leaving a total of 17 globular clusters. Our final sample con- 
tained the following clusters: NGC 1904, 2298, 5024, 5053, 5466, 5897, 6101, 
6205, 6254, 6341, 6397, 6535, 6809, 7078, 7099, 7492, and Terzan 8. 

One of the things we checked is the the inferred dispersion in the age 
of the 17 globular clusters is not larger than that which we expect based 
on the observational uncertainties. Using the uncertainties in the individual 
ages determined from uncertainties in the observed turnoff magnitude and 
metallicity, we examined the dispersion about the mean age for the 17 clusters 
using a test. We found a reduced of 0.55 per degree of freedom, 
indicating both no evidence for any intrinsic dispersion in age for our sample, 
and that the quoted observational uncertainties for each cluster may be too 
generous. In any case, given the quoted accuracy, it is certainly consistent 
to assign a single mean age for the sample. 

One chief observational uncertainty common to all globular clusters is 
retained until the end of the analysis. This is the determination for 
RR-Lyrae variables. In order to determine AV^b as a function of age and 
metallicity, we combine our theoretical turn-off magnitude with an observa- 
tionally based estimate for the absolute magnitude of the HB, in the RR Lyr 
instability strip (hereafter referred to as Mv(RR)). There are a number of 
independent, observationally based techniques which can be used to derive 
Mv(RR). In general, it has been found that the absolute magnitude of the 
RR Lyr stars can be represented by an equation of the form 

Mv(RR) = /i[Fe/H]+7, (3) 

where /i is the slope with metallicity and 7 is the zero-point. Note that 
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Mv(RR) is independent of age (at least, with systems greater than 8 Gyr 
old). Recent estimates for the slope with metallicity vary from 0.15 to 0.30 
pB[ . Fortunately, since we are determining the mean age of 17 globular 
clusters in the restricted metallicity range —2.41 < [Fe/H] < —1.60, the 
uncertainties in the slope has only a small effect on our age estimate. Tests 
which we conducted indicated that the maximum difference in our mean age 
was only 0.5% when the slope was varied between 0.15 and 0.30. The most 



recent work suggests that /i = 0.20 is likely to be correct ^ which is the 
value we have used in this study. Uncertainties in the zero-point, 7 in eq. ^ 
have a large impact on our derived age estimates. For this reason, we have 
spent considerable time reviewing recent observational estimates of the zero- 
point. These estimates for the zero-point are usually given as a Mv(RR) value 
at a specific metallicity. As the globular clusters in our sample have a median 
metallicity of [Fe/H] = —1.82 and a mean metallicity of [Fe/H] = —1.93 we 
have elected to transform the various zero-point estimates to [Fe/H] = —1.90 
using a slope oi ^ = 0.20 ± 0.04. 

Layden Hanson & Hawley [^ have recently used the statistical parallax 
technique to determine Mv(RR) = 0.68 ± 0.12 in field halo RR Lyr stars. 
Walker [^ measured the apparent magnitude of LMC RR Lyr stars, and 
assumed an LMC distance modulus of 18.5 to infer Mv(RR) = 0.44 ± 0.10. 
This choice for the distance modulus of the LMC was based on the Cepheid 
distance, main sequence fitting and the SN1987A ring distance to the LMC. 
This last method is a purely geometrical method, and should be the most 
reliable. However, the SN1987A distance to the LMC has recently been re- 
vised to 18.37 [0, implying Mv(RR) = 0.57 ± 0.10. Main sequence fitting 
of globular cluster CM diagrams to local halo stars with well determined 
parallaxes can be used to determine the distance to globular clusters, and 
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hence, Mv(RR). Unfortunately there is only one relatively metal-rich sub- 
dwarf which has a well determined parallax. Application of this technique to 
the globular cluster M5 yields M^(RR) = 0.76 ±0.12 0. The only direct 
determination of Mv(RR) in a metal-poor globular cluster is by Storm, Car- 
ney & Latham |Q. They used a Baade-Wesselink/infrared flux analysis to 
determine Mv(RR) = 0.52 ±0.26. Although the error is large due to possible 
systematic uncertainties, it does suggest that the RR Lyr stars in metal-poor 



globular clusters are somewhat brighter than those found in the field , or 



in metal- rich globular clusters ||35|. In light of the above estimates, we have 
elected to use Mv(RR) = 0.60±0.08 (corresponding to 7 = 0.98±0.08). This 
central value was chosen by a straight average of the 4 published Mv(RR) es- 
timates referenced above. The error bar was chosen to ensure that the 1 a 
range would include the central value obtained for Mv(RR) in a metal-poor 
globular cluster, and the 2 a range (0.44 - 0.76) would encompass all of the 
estimates quoted above. Our choice is further supported by a study which 
just appeared comparing four different distance estimators, including kine- 
matic, RR Lyr, Cepheid and Type II Supernovae for consistency [^, and 
found that this appeared to require a range for Mv(RR) similar to the one 
we chose. 

Figure 1 displays the ensemble of age estimates from our Monte Carlo, 
for different values of Mv(RR). It was generated as follows: For each of the 
sets of isochrones and given a value for Mv(RR) we determine a mean age 
and \a uncertainty in the mean. In the figure we show these values for each 
isochrone set, for 3 values of Mv(RR): the mean, .60, and the endpoints of 
our 95% C.L. range, .44 and .76. In order to obtain the final histograms 
displayed as figs 2 a-b, we follow an analogous procedure, but this time allow 
Mv(RR) to be a random variable, and sample the sets of isochrones with 
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replacement 12,000 times. For each sample, rather than the mean age, we 
record a random age drawn from a gaussian distribution with the mean age 
and variance for that isochrone set at that Mv(RR). Then the data are sorted 
and binned to produce the figures: 2(a) the full distribution for the assumed 
Gaussian spread in Mv(RR) of ±.08, and 2(b) the distribution for fixed value 
of Mv(RR) of 0.6. In this way the effect of the uncertainty in Mv(RR) can 
be explicitly examined. 

Conclusions 

Our results indicate that at the one sided 95% confidence level (deter- 
mined by requiring 95% of the determined ages to fall above this value) a 
lower limit of approximately 12.1 Gyr can be placed on the mean value of 
these 18 Globular clusters. (The symmetric 95% range of ages about the 
mean value of 14.56 Gyr is 11.6-18.1 Gyr.) Note that the distribution de- 
viates somewhat from Gaussian, as one might expect. In particular, at the 
lower age limits the rise is steeper than Gaussian, reflecting the fact that 
essentially all models give an age in excess of 10 Gyr, while the tail for larger 
ages is larger than gaussian. The explicit effect of the largest single common 
observational uncertainty, that in Mv(RR), increases the net width of the 
distribution by approximately ±0.6 Gyr (i.e. f« ±5%). Note that simply 
varying Mv(RR)over its full 2a range, keeping all other parameters fixed, 
would produce a ±16% change in GG ages estimates. For comparison, the 
next most significant input parameter uncertainties in this same sense are 
[a/Fe] (±7% effect), mixing length (±5% effect), and diffusion, -"-^Np reac- 
tion rate, and primordial Helium abundance, each of which would affect age 
estimates at the ±3% level if allowed to vary over its entire range, keeping 
all other parameters fixed. 
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We believe that this result can now be used with some confidence to 
compare to cosmological age estimates. Of course, in addition to the age 
determined here one must add some estimate for the time it took our galactic 
stellar halo to form from the initial density perturbations present during the 
Big Bang expansion. Estimates for this formation time vary from 0.1 - 
2 Gyr. To be conservative, we choose the lower value. In this case, we 
find that the age of globular clusters in our galaxy is inconsistent with a 
fiat, matter dominated universe unless h < 0.54, and for a nearly empty, 
matter dominated universe unless h < 0.80. If the value of h is definitely 
determined to be larger than cither of these values, some modification, such 
as the addition of a cosmological constant, would seem to be required. 
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Figure Captions 

Figure 1: The suite of Models generated by the Monte Carlo Procedure. For 
each Model the 1 — a age range is plotted for 3 choices of Mv(RR). The red, 
green and blue data points show the age variation for Mv(RR) = 0.44, 0.60 
and 0.76 respectively. 

Figure 2: The histograms exhibit the relative numbers of realizations of mean 
globular cluster ages drawn randomly from the Monte Carlo data set (with 
uncertainties on individual age estimates taken to be Gaussian) using a value 
for Mv(RR) (the absolute magnitude of the RR Lyrae Variables) chosen from 
one of two different distributions: a) a Gaussian, Mv(RR) = 0.6 ± .08, b) a 
Delta function, Mv(RR) = 0.6. The dashed line is a Gaussian approximation 
to the actual distribution. The mean and standard deviation of the Gaussians 
are shown in the Legend. The horizontal arrows show the 1 and 2-a ranges in 
age based upon the Gaussian approximation. The one-sided 95% Confidence 
Limit for a lower bound on the age of the Universe is displayed as an arrow 
extending to the right from a vertical bar. This is calculated directly from 
the generated distribution. 
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